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We propose methods for constructing high-quality pseudorandom number generators (RNGs) 
based on an ensemble of hyperbolic automorphisms of the unit two-dimensional torus (Sinai-Arnold 
map or cat map) while keeping a part of the information hidden. The single cat map provides the 
random properties expected from a good RNG and is hence an appropriate building block for an 
RNG, although unnecessary correlations are always present in practice. We show that introducing 
hidden variables and introducing rotation in the RNG output, accompanied with the proper initial- 
ization, dramatically suppress these correlations. We analyze the mechanisms of the single-cat-map 
correlations analytically and show how to diminish them. We generalize the Percival- Vivaldi theory 
in the case of the ensemble of maps, find the period of the proposed RNG analytically, and also 
analyze its properties. We present efficient practical realizations for the RNGs and check our pre- 
dictions numerically. We also test our RNGs using the known stringent batteries of statistical tests 
and find that the statistical properties of our best generators are not worse than those of other best 
modern generators. 

PACS numbers: 02.50.Ng, 02.70.Uu, 05.45.-a 



I. INTRODUCTION 

Molecular dynamics and Monte Carlo simulations are important computational techniques in many areas of science: 
in quantum physics 1] , statistical physics 2] , nuclear physics "s^ , quantum chemistry 4] , material science , among 
many others. The simulations rely heavily on the use of random numbers, which are generated by deterministic 
recursive rules. Such rules produce pseudorandom numbers, and it is a great challenge to design random number 
generators (RNGs) that behave as realizations of independent uniformly distributed random variables and approximate 
"true randomness" 

There are several requirements for a good RNG and its implementation in a subroutine library. Among them are 
statistical robustness (uniform distribution of values at the output with no apparent correlations), unpredictability, 
long period, efficiency, theoretical support (precise prediction of the important properties), portability and others [a, 
it 

A number of RNGs introduced in the last five decades fulfill most of the requirements and are successfully used in 
simulations. Nevertheless, each of them has some weak properties which may (or may not) influence the results. 

The most widely used RNGs can be divided into two classes. The first class is represented by the Linear Gongruential 
Generator (LCG), and the second, by Shift Register (SR) generator. 

Linear Congruential Generators (LGGs) are the best-known and (still) most widely available RNGs in use today. An 
example of the realization of an LGG generator is the UNIX rand generator ?/„ = (1103515245 j/n-i + 12345) (mod 
2'^-'^). The practical recommendation is that LCGs should be avoided for applications dealing with the geometric 
behavior of random vectors in high dimensions because of the bad geometric structure of the vectors that they 
produce 

Generalized Feedback Shift Register (GFSR) sequences are widely used in many areas of computational and simula- 
tional physics. These RNGs are quite fast and possess huge periods given a proper choice of the underlying primitive 
trinomials p!o|. This makes them particularly well suited for applications that require many pseudorandom numbers. 
But several flaws have been observed in the statistical properties of these generators, which can result in systematic 
errors in Monte Carlo simulations. Typical examples include the Wolff single cluster algorithm for the 2D Ising 
model simulation random and self-avoiding walks and the 3D Blume-Capel model using local Metropolis 
updating [l^ . 

Modern modifications and generalizations to the LCG and GFSR methods have much better periodic and statistical 
properties. Some examples are the Mersenne twister |3| (this generator employs the modified and generalized GFSR 
scheme) , combined LCGs generators [iBj and combined Tausworthe generators 0, 0| . 

Most RNGs used today can be easily deciphered. Perhaps the generator with the best unpredictability properties 
known today is the BBS generator Il9j|. which is proved to be polynomial-time perfect under certain reasonable 
assumptions i 0| if the size s of the generator is sufficiently large. This generator is rather slow for practical use 
because its speed decreases rapidly as s increases. The discussion of cryptographic RNG is beyond our analysis. 
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We propose using an ensemble of simple nonlinear dynamical systems to construct an RNG. Of course, not all 
dynamical systems are useful. For instance, baker's transformation is a simple example of a chaotic system: it is area 
preserving and deterministic, and its state is maintained in a bounded domain. The base of baker's transformation is 
the Bernoulli shift Xn+i = 2x„ ( mod 1), it yields a sequence of random numbers provided we have a random irrational 
seed. But in real computation, the seed number has finite complexity, and the number of available bits decreases at 
each step. Obviously, there is no practical use of this scheme for an RNG. 

The logistic map [23, E| also does not help to construct an RNG. First, manipulation with real values of fixed 
accuracy leads to significant errors during long orbits. Second, the sequence of numbers generated by a logistic map 
does not have a uniform distribution |2j. Also, the logistic map represents a chaotic dynamical system only for 
isolated values of a parameter. Even small deviations from these isolated values lead to creating subregions in the 
phase space, i.e., the orbit of the point does not span the whole phase space. 

The next class of dynamical systems is Anosov diffeomorphisms of the two-dimensional torus, which have attracted 
much attention in the context of ergodic theory. Anosov systems have the following stochastic properties: ergodic- 
ity, mixing, sensitive dependence on initial conditions (which follows from the positivity of the Lyapunov exponent), 
and local divergence of all trajectories (which follows from the positivity of the Kolmogorov-Sinai entropy). These 
properties resemble certain properties of randomness. Every Anosov diffeomorphism of the torus is topologically 
conjugate to a hyperbolic automorphism, which can be viewed as a completely chaotic Hamiltonian dynamical sys- 
tem. Hyperbolic automorphisms are represented by 2x2-matrixes with integer entries, a unit determinant, and real 
eigenvalues, and are known as cat maps (there are two reasons for this terminology: first, CAT is an acronym for Con- 
tinuous Automorphism of the Torus; second, the chaotic behavior of these maps is traditionally described by showing 
the result of their action on the face of the cat H^l). We note that cat maps are Hamiltonian systems. Indeed, if 
k = Tr(M) = mil +^22 > 2, then the action of map on the vector can be described as the motion in the phase 

space specified by the Hamiltonian |2^ H{p, q) — (fc^ — 4)"^/^ sinh^^((fc^ — 4)^/^/2)(r7ii2p^ — m2iq'^ + (mn — TO22)m)- 
Here, p and q are taken modulo 1 at each observation (i.e., we preserve only the fractional part of p and q; the integer 
part is ignored), and observations occur at integer points of time. 

In this paper, we present RNGs based on an ensemble of cat maps and analyze the requirements for a good RNG with 
respect to our scheme. The basic idea is to apply the cat map to a discrete set of points (there are two modifications: 
for (7 = 2™ and for prime g, where (gxg) is the lattice) such that each point belongs to a different periodic trajectory. 

A similar utilization of cat maps for an RNG is called the matrix generator for pseudorandom numbers. It was 
introduced in j24li.l25j | and discussed for prime values of g. But because the single matrix generator is a generalization 
of the Hnear congruential method, it suffers from both the defects of LCG |2^ and the defects of GFSR (see Sec. lIVp . 
The periodic and statistical properties of the matrix generators and of the equivalent multiple recursive generators 
have been studied [27^ .28,] , but the single 2 x 2-matrix generator still has significant correlations between values at the 
output. 

Also, there is an impressive theoretical basis for relating properties of the periodic orbits of cat maps and properties 
of algebraic numbers j29l | , which to the best of our knowledge has never been directly applied to RNG theory. Applying 
the ensemble of matrix transformations of the two-dimensional torus while using only a single bit from the point of 
each map, and utilizing rotation in the RNG output are the distinctive features of our generator. Also, as for other 
generators, a proper initialization of the initial state is important. As will be seen, the proposed scheme has several 
advantages. First, it can essentially reduce correlations and lead to creating an RNG not worse than other modern 
RNGs. Second, both the properties of periodic orbits and the statistical properties of such a generator can be 
analyzed both theoretically and empirically. Several examples of RNGs made by this method, as well as the effective 
realizations, are presented. 

The generator is introduced in Sec.m In Sec. lIIII we present the results for stringent statistical tests. Correlations 
for a single cat map are also analyzed thoroughly, and some correlations are found by the random walks test. We 
analyze the mechanism of these correlations in Sec. II VI they appear to be associated with the geometric properties of 
the cat map. We find these correlations analytically fSec. IIV|I . We provide a method for obtaining quantities such as 
the periods of cat maps, the number of orbits with a given period, and the area in the phase space swept by the orbits 
with a given period (Appendix ^ . We also provide a method for obtaining periods of the generator for arbitrary 
parameters of the map and lattice f Appendix IB|) . This gives the primary theoretical support of the generator. In 
particular, we find that the typical period of the generator for the 2™x2™ lattice is = 3 • 2™"^. The method is 
based on the work of Percival and Vivaldi, who transformed the study of the periodic orbits of cat maps into the 
modular arithmetic in domains of quadratic integers. The key ideas needed for our consideration are briefly reviewed in 
Appendix^ Appendix [C] gives the method for analyzing correlations between orbits of different points and choosing 
the proper initial conditions to minimize the correlations. AppendixiDland supplemental details ^45j support the other 
sections, giving detailed proofs of the underlying results. Appendix |E| presents the efficient realizations for several 
versions of RNG, the initialization techniques, and the analysis of the speed of the RNGs. 
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II. THE GENERATOR 



A. Description of the method 



We consider hyperbolic automorphisms of the unit two-dimensional torus (the square (0, 1] x (0, 1] with the opposite 
sides identified). The action of a given cat map R is defined as follows: first, we transform the phase space by the 
matrix 



Af= "^'A eSMZ); (f) 

\^ TO21 m,22 J 

second, we take the fractional parts in (0, f) of both coordinates. Here S'L2(Z) denotes the special linear group of 
degree 2 over the ring of integers, i.e., the elements of M are integers, detM = 1, and the eigenvalues of M are 
A = (fc± Vfc^ — 4)/2, where k = Tr(M) is the trace of the matrix M. The eigenvalues should be real because complex 
values of A lead to a nonergodic dynamical process, and the hyperbolicity condition is > 2. 

It is easy to prove that the periodic orbits of the hyperbolic toral automorphism R consist precisely of those points 
that have rational coordinates [S^, I23L |29| . Hence, it is natural to consider the dynamics of the map defined on the 
set of points with rational coordinates that share a given denominator g. The lattice of such points is invariant under 
the action of the cat maps. In practice, we construct generators with g — 2™, where m is a positive integer, and 
generators with g = p = 2™ — 1, where m is a Mersenne exponent, i.e., p = 2™ — 1 is a prime. 

The notion of an RNG can be formalized as follows: a generator is a structure G — {S, so,T,U,G), where 5' is a 
finite set of states, sq € 5* is the initial state (or seed), the map T : S S is the transition function, U is a finite 
set of output symbols, and G : 5 ^ C/ is the output funet^on U ■ 

Thus, the state of the generator is initially sq, 
and the generator changes its state at each step, calculating s„ = T(s„_i), u„ = G(s„) at step n. The values u„ 
at the output of the generator are called the observations or the random numbers produced by the generator. The 
output function G may use only a small part of the state information to calculate the random number, the majority 
of the information being ignored. In this case, there exist hidden variables, i.e., some part of the state information is 
"hidden" and cannot be restored using only the sequence of RNG observations. 

We consider the generator with S* = L®, where L = {0, 1, . . . , 5 — 1} x {0, 1, . . . , 5 — 1} is the lattice on the torus 
and s is a positive integer. In other words, the state consists of coordinates of s points of the gxg lattice on the torus. 

/x^^^^\ fni (0) 

For instance, the initial state consists of points ( (o)) , where x\ ,yl £ {0, 1, . . . , g — 1} and i = 0, 1, . . . , (s — 1). We 

note that these are points of the integer lattice, i.e., xf^ and yf^ are positive integers. The actual initial points on 
the unit two-dimensional torus (0, 1] x (0, 1] are 



4^19 

yf'/9. 



i = 0,l,...,(s-l). (2) 



The transition function of the generator is defined by the action of the cat map R, i.e., these s points are affected 
at every step by the cat map: 

*=o,i,...,(.s-i). (3) 

Vi 19/ KVi 7.9/ 

Here the mod 1 operation means taking the fractional part in (0, 1) of the real number. An equivalent description of 
the transition function is 

x) \ fx: 



(«) in-i) ) (mod g), z = 0, 1, . . . , (s - 1). (4) 

We let a^^'' denote or 1 depending on whether x^-^-' < {g/2) or > {g/2), i.e., a I") = [2x,^"VgJ. The output 
function of the generator G : ^ {0, 1, . . . , 2^ — 1} is defined as a^"^ = J^iZo Q^!"'' ' 2^ In other words, a^"^ is an 
s-bit integer consisting of the bits ao""*, a^""*, . . . , a["^i. In the case g — 2"^, a*^"^ contains precisely the first bits of the 
integers Xq"^ . . . , x^J^^^. The sequence of random numbers produced by the generator is {a^"''}. 
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We see that the constructed RNG has much hidden information. For example, ii g = 2™, then s{m~l) bits of (^(„)) 

are the hidden variables; these are the bits that are not involved in constructing the value of the output function a*^"^ 
Thus, applying the chaotic behavior of Anosov motion and introducing an ensemble of systems while keeping part of 
the information hidden are the main ingredients of the proposed method. Good stochastic properties of the underlying 
continuous system are obviously necessary for good generators. For example, the logarithm of the multiplier in the 
continuous transformation of the LOG can be viewed as the Lyapunov exponent, which is always greater than 1, and 
this leads to the divergence of trajectories. The huge number of points on a lattice makes the continuous system a 
good first approximation to the RNG and leads to the importance of good chaotic properties. Introducing hidden 
variables reduces correlations (as is shown in Sec. 

The calculation of the period of the RNG is presented in Appendix IbI The typical period length is T„i = 3 • 2™^^ 
for the 2'"x2'" lattice. The proper initializations for the generators are presented in Appendix IeI The proper 

initialization guarantees that the actual period is not smaller than Tm and that the points ( (o)), i = 0, 1, . . . , (s — 1), 

belong to different orbits of the cat map. 



B. Connection with other generators 

There are several known connections between Anosov dynamical systems and pseudorandom number generation. 

First, the concept of the Shift Register Sequence, which is widely used to construct high-quality RNGs, is con- 
nected to dynamical systems (see, e.g., the discussion in |2l|)- Let the state of the shift register be Vn-i = 
(a„_r, a„_r+i, • • ■ , flri-i)- At the next iteration, the state of the shift register is Vn — (a„_r+i, an-r+2, • ■ • , an), 
where a„ — Cr-an-r + Cgfln-s (mod 2). In other words, Vn+i = Av^ (mod 2), where A is an (rxr)-matrix. 

Second, LCGs in some cases can be described by the action of the hyperbolic toral automorphism [37l |. 

Last, it can be shown that, for each i, the sequence {x^"''}, defined above, as well as the sequence {y|"''}, follows a 
linear recurrence modulo g: 



= fcx^"-!) - ga;("-2) (mod g) (5) 
yin) ^ ky{n-l) _ qy("-2) (inod g), (6) 

where fc = Tr (M) , and q = det M = 1. The characteristic polynomial of the last linear recurrence is f{x) = q—kx+x^, 
which is exactly the same as that of the matrix M ||2J, . 

The period properties of sequence jSJ follow from the arithmetical methods for q — 1 (see Appendix ^ and 
Appendix ^ and from the finite field theory in the case where g = p is a prime . 



C. Generators for prime g: modifications for detM = 1 and for det Af 7^ 1 

The matrix generator of pseudorandom numbers equivalent to sequence Q was studied in |^ 0, Ull in the case 
where g = det M = p is a prime. Sequence lO yields the maximum possible period if and only if the characteristic 
polynomial f{x) is primitive over Zp. But for q — 1 the polynomial f{x) = x"^ — kx + 1 is not primitive over Zj, for 
p > 2. Therefore, if detM = 1, the period is always smaller than — 1. An even stronger result follows from (29j: 
the period cannot be larger than p + 1 when g is a prime and det M — I. 

Matrix generators with q 1 are not immediately connected with Hamiltonian dynamical systems. Indeed, the 
transformation with det M ^ 1 does not preserve the volume in phase space and does not immediately represent a cat 
map. However, whatever q is, we have detM^"^ = 1 (mod p). This means that the action of the matrix M'p~^ on a 
lattice p X p is exactly the same as the action of a unimodular matrix. Therefore, any orbit of a "non-Hamiltonian" 
transformation M contains exactly p — 1 cat-map orbits. 

Also, transformations with q — 1 preserve the norm on the orbit modulo g (see Appendix Q, in contrast to 
transformations with q^ \. It is shown in Appendix O that some of the correlations between the orbits are inherent 
in the case q = 1 and are suppressed for q ^ 1. 



D. Rotating the RNG output 

It will be seen that in the scheme of the generator, there are correlations between the first bits of 0^"^ correlations 
between the second bits of 0*^"^ and so on. To suppress these correlations, we modify the algorithm as follows. At each 
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step, we renumber the points in the generator output: 1 ^ 2, 2 — > 3, . . . , s — > 1. In other words, the bits inside a*^") are 
rotated, and the RNG output function is defined as = Y.tZo ■ 2^*+") ^ """'^ instead of a^") = Y.tZo a*"'* ' 2% 
where aj-""* = [2x|"-'/(7j . 

The main advantage of the modified algorithm is that it leads to decreasing the correlations of the values a'"' 
between each other. For example, we will see in Sec. IIVI that the rotation strongly reduces the specific correlations 
found by the random walks test. 

We note that the rotating the bits in the RNG output does not deteriorate any properties of the RNG provided 
that s divides the period of free orbits Tm (in practice, this is a very realistic condition). In particular, neither does 
the generator period become smaller (see Appendix IbI. nor do the statistical properties become worse. 

Rotating the bits in the RNG output is thus a practically useful modification. In addition, the rotation makes 
deciphering an even more complicated problem. 

III. STATISTICAL TESTS 
A. Simple Knuth tests 

In this section, we present the results of several standard statistical tests 6] that reveal the correlation properties 
of the generator described in Sec. Namely, the frequency test, serial test, maximum-of-t test, test for monotonic 
subsequences ("run test") and collision test were applied for an RNG with M = (g g), 5 = 2™ = 2^* and s = 28 
points in the state. All the statistical tests were passed. All empirical tests except the collision test (CT) are based 
on either the chi-square test (x^) or the Kolmogorov-Smirnov test (KS). We follow Knuth's notation d\. 

The results of the tests are presented in Table U where n is the number of values of a'"^ for each test (for the serial 
test, n is the number of pairs {a'^"-', a*^^"+"'^^}) and v is the number of degrees of freedom. For the serial test d — &, 
i.e., we used exactly 3 bits of each a*-"-* number; hence, v — cP ~ 1 — 63. For the run test, v — 5 means that we sought 
monotonic subsequences of lengths 1,2,3,4,5 and of length > 6. 

For all of the KS tests, the empirical distributions of P(A'+) and P{K~) were calculated, where P{x) is the 
theoretical Kolmogorov-Smirnov distribution Figure ^ shows these empirical distributions for the frequency test. 
These distributions lead to their own values of and K~: the values P{K''^) and P{K'~) characterizing the 
empirical distribution of and the values P{K'^) and P{K'~) characterizing the empirical distribution of K~ . 
These values are presented in TableQ] Our RNG passes all the KS tests because the values Ar+ and are distributed 
in accordance with the theory prediction. For each chi-square test, the empirical distribution of 20 values of P{V) was 
calculated. In our tests, it looks similar to those shown in Fig.^ where P{x) is the theoretical chi-square distribution 
and V is the output of the chi-square test. For each collision test, the number of collisions c and the theoretical 
probability P(c) that the number of collisions is not larger than c were calculated. The empirical distribution of P{c) 
was analyzed, and the results are also presented in Table 



TABLE L Results of the statistical tests for the RNG based on the ensemble of cat maps (see Sec^J with parameters M = (^ g) , 
3 = 2™ = s = 28. 



Test 


Parameters 


Number and 
type of tests 


Tests output values 
Vi and V2 


Distribution of Vi 


Distribution of V2 


Conclusion 


P{K'+) 


P{K'-) 


P{K'+) 


P{K'-) 


Frequency test 


n= 10® 


20 


KS 


Vi^K+,V2 = K- 


0.592392 


0.174451 


0.457758 


0.473379 


PASSED 


Serial test 


n = 10^(^ = 8 


20 




Vx = V 


0.837112 


0.17128 


n/a 


n/a 


PASSED 


Run test 


n = l(f,u^ 5 


20 




Vi = V 


0.383601 


0.805434 


n/a 


n/a 


PASSED 


Maximum-of-t test 


n= 10^^ = 5 


20 


KS 


Vi =K+,V2 = K- 


0.120912 


0.765201 


0.704026 


0.589702 


PASSED 


Collision test 


m = 22°,n = 2" 


20 


CT 


Vi=c 


0.150537 


0.858955 


n/a 


n/a 


PASSED 



We note that all the empirical tests here except the collision test are essentially multibit. This means that the 
whole ensemble of cat maps influences the test result, and one can guess that hidden variables inside the generator 
is one reason for the successful test results. A single-bit cat-map generator, i.e., a generator from Sec. ^ with s = 1, 
does not contain hidden variables. Most of the tests for a single-bit cat-map generator are also successfully passed. 
Namely, the frequency test and the serial test, which were modified for a one-bit generator, and the collision test are 
passed. But there are correlations in the single-bit cat-map generator (discussed later in this paper), and the most 
convenient method for observing them is the random walks test with = 1/2 (see Sec. IIV|) . The random walks test 
is not the only test that can reveal the single-bit cat-map correlations. The same correlations are also observed by 
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K* 



0.0 0.2 0.4 0.6 0.8 1.0 

K- 



FIG. 1: Distribution of P{K^) and P{K ) for the frequency test. The test was performed for the RNG from Sec.|n]with the 
parameters M ^ Ql) , g ^ 2"^ = 2'^^ , s = 28. 

improved versions of some of the standard tests, e.g., the serial test for subsequences of length 5. Of course, many 
tests in Sec. IIIIBl would not be passed by a single-bit cat-map generator. 

For comparison, we analyzed a simple generator based on the single cat map. Tabled shows that such a generator 

with the transition function defined as (^(„)) = M(^(„_i)) (mod 1) and the output function defined ^"^ has 

very bad properties. Of course, the frequency test is passed, since the trajectories of a cat map uniformly fill the 
phase space. But all the other tests are failed. Therefore, the simple generator based on the single cat map does have 
strong correlations in the output and is not useful practically. 



TABLE IL Results of the statistical tests for a simple RNG based on the single cat map with parameters M = (^ |?) , g = 2"" = 



Test 


Parameters 


Number and 
type of tests 


Tests output 
values Vi, V2 


Distribution of V\ 


Distribution of V2 


Conclusion 


P(K'+) 


P{K'~) 


P{K'+) 


P{K'-) 


Frequency test 


n = 10'' 


20 


KS 


K+;K- 


0.96235 


0.170581 


0.787296 


0.067341 


PASSED 


Serial test 


n = 5- 10^d = 8 


20 




V 


0.989499 


0.006013 


n/a 


n/a 


FAILED 


Run test 


n = 10^^/ = 5 


20 


x' 


V 





1 


n/a 


n/a 


FAILED 


Maximum-of-t test 


n = 10^^ = 5 


20 


KS 


K+;K- 





1 





1 


FAILED 


Collision test 


d = 4,m = 2"',n = 2^'' 


20 


CT 


c 





1 


n/a 


n/a 


FAILED 



B. Batteries of stringent statistical tests 

Knuth tests are very important but still not sufficient for the present-day sound analysis of the RNG statistical 
properties. Hundreds of statistical tests and algorithms are available in software packages, for example, widely used 
packages DieHard |3^J, NIST [s^ and TestUOl j34j. All of them include tests, described by Knuth Q, as well as many 
other tests. 

Table IIIII shows the summary results for the SmallCrush, PscudoDiehard, Crush and Bigcrush batteries of tests 
from 34]. SmallCrush, PseudoDiehard, Crush and Bigcrush contain 14, 126, 93 and 65 tests respectively. The detailed 
parameters and initializations for the generators GS, GR, GSI, GRI, GM19 and GM31, based on the scheme proposed 
in Sec.^ are given in Appendix IeI 

For comparison, we also test several other generators, namely, the standard generators RAND, RAND48 and 
RANDOM and the modern generators MT19937, MRG32k3a and LFSR113. RAND is the simple LCG generator 
based on the recursion Xn = (1103515245 x„_i + 12345) (mod 2^^). RAND48 is the 64-bit LCG based on the 
recursion x,-, = 25214903917 a;„_i + 11 (mod 2'**). RANDOM provides an interface to a set of five additive feedback 
random number generators. RAND, RAND48 and RANDOM are implemented in the functions rsindO, rand48() 
and randoniO in the standard Unix or Linux C library stdlib (see the documentation to randO, rand48() and 
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randomO). MT19937 is the 2002 version of the Mersenne Twister generator of Matsumoto and Nishimura 0|, which 
is based on the recent generahzations to the GFSR method. MRG32k3a is the combined multiply recursive generator 
proposed in and LFSR113 is a combined Tausworthe generator of L'Ecuyer [T^ . 

The detailed statistics for the batteries of tests and the explicit results for every single test from the batteries can 
be found in |4^ . 

TABLE III: Numbers of failed tests for the batteries of tests SmallCrush, Crush, Bigcrush |3^ . and DieHard [3fll |. Here 
k = Tr(M) and q — detM are the RNG parameters (see Sec. UTI and Appendix IeI. For each test, we present three numbers: 
the number of statistical tests with p-values outside the interval [10~^, 1 — 10^^], number of tests with p-values outside the 
interval [10~^, 1 — 10^^], and number of tests with p-values outside the interval [10"^'', 1 — 10~^°]. 



Generator 


k 


q 


SmallCrush 


Diehard 


Crush 


Bigcrush 


GS 


3 


1 


0,0,0 


44, 29, 29 


20, 16, 14 


22,20, 19 


GR 


3 


1 


0,0,0 


5,0,0 


5,1,0 


15,10,7 


GSI 


11 


1 


0,0,0 


1,0,0 


10,1,0 


13, 7,6 


GRI 


11 


1 


1,0,0 


6,0,0 


5,0,0 


13,6,5 


GM19 


15 


28 


0,0,0 


2,0,0 


2,0,0 


3,0,0 


GM31 


7 


11 


0,0,0 


2,0,0 


3,0,0 


1,0,0 


RAND 






13,13,12 


88, 84, 82 


102, 100, 100 


85, 83, 79 


RAND48 






5,5,3 


27, 23, 22 


22,20,20 


27, 23, 22 


RANDOM 






3,2,2 


17,15,15 


13,11,10 


21, 15, 14 


MRG32k3a 






1,0,0 


3,0,0 


4,0,0 


2,0,0 


LFSR113 






0,0,0 


3,0,0 


8,6,6 


8,3,3 


MT19937 






0,0,0 


2,0,0 


1,0,0 


4,0,0 



We consider the test "failed" if the p- value lies outside the region [10~^, 1 — 10~^]. Most of the p-values for the 
failed tests for the cat-map generators are of the order of 10"'^ to 10~^, but several are very small. We believe that the 
reason for small p-values is connected with the small period of the generators GS, GR, GSI, GRI and UNIX RAND. 
A period of the order of 3- 10^, while sufficient for some applications, is not sufficient for many of the tests from Crush 
and Bigcrush. Therefore, the generators GS, GR, GSI and GRI demonstrate smaller p-values and larger numbers of 
failed tests from Crush and Bigcrush. 

The existence of linear congruential dependences between orbits is another reason for small p-values for GS, GR, 
GSI and GRI. These correlations are described analytically in Appendix El The GM19 and GM31 generators, 
having a period sufficient for the Crush and Bigcrush batteries, are simultaneously free from the linear congruential 
dependences. Therefore, they demonstrate much better statistical properties in Table ITTll 

Because we apply hundreds of tests, the number of failed tests is susceptible to random statistical flukes, especially 
when the p-values of failed tests lie in the suspect region [10~^, 10~^] U [1 — 10~^, 1 — 10~^]. Table HVl illustrates the 
flukes by showing the results of all batteries of tests for the generator GM31. The batteries were executed in the order 
SmallCrush, SmallCrush, PseudoDiehard, PseudoDiehard, Crush, Crush, BigCrush and BigCrush, i.e. each battery 
was executed twice. For the tests in Tables Hill and Hvl the generator GM31 was initialized with identical parameters, 
in accordance with Appendix El The numbers of failed tests themselves in Table IIIII and Table IIVI approximately 
indicate the statistical robustness of the generators. But if the p-value lies in the suspect region, one does not know 
exactly whether systematic correlations were found in the RNG or a statistical fluke occured. 

We conclude that the best of the generators based on cat maps are competitive with other good modern generators. 
In particular, we recommend the RNG realizations for GRI and GM31 for practical use. In Appendix |e1 we present 
the effective realizations of the generators GRTSSE and GM31-SSE and recipes for the proper initialization. Among 
the generators examined here these are the best respective realizations with g — 2™ and with prime g. 

IV. THE RANDOM WALKS TEST 

Analyzing the statistical properties of the generator theoreticall y is another important challenge. Such an analysis 
is traditionally performed by discussing the lattice structure 12^1281 Isof and discussing the discrepancy [s^l . The 
discrepancy of a matrix generator was analyzed by Niederreiter |27l |. who in particular, proved that the behavior of 
the discrepancy is strongly connected to the behavior of an integer called the figure of merit. Although calculating 
the exact values of the figure of merit would give an excellent basis for the practical selection of matrices for matrix 
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TABLE IV: Applying each battery of tests twice for tlie GM31 generator. Here the test is considered failed if the p-value lies 
outside the interval [10"^, 1 - 10"^] 





WiiTTiriPT' nf 


Failed tests, testing first time 


Failed tests, testing second time 




icLlltLl ±fc!bLb 




N 3,1116 


p- value 


IN O 


Name 


p- value 


oiiiaiiv^i Lihii 


n /n 














PseudoDiehard 


o/ Z 


1 


BirthdaySpacings 




r- 
D 


CoUisionOver 


U.UUi4 






D 


CollisionOvGr 


n QQQ 


7 


CollisionOver 


n nm s 

U.UUio 






14 


Run of UOl 


0.0093 








Crush 


2/2 


70 


Fourierl, r = 


0.0095 


14 


BirthdaySpacings, t — 7 


0.0024 






71 


Fourier 1, r = 20 


0.0093 


70 


Fourierl, r = 


0.0033 


BigCrush 


4/3 


6 


MultinomialBitsOver 


0.0014 


32 


SumCoUector 


0.0067 






23 


Gap, r — 


0.9970 


36 


RandomWalkl J(L=90) 


0.9952 






27 


CoUisionPermut 


0.0052 


41 


Random Walkl J(L=10000) 


0.9965 






39 


Random Walkl H(L=1000) 


0.9982 









generators, these values are still very hard to compute. To the best of our knowledge, this calculation has never been 
done for matrix generators of pseudorandom numbers. 

Because of the hidden variables, the lattice structure of the matrix generator does not directly influence the statistical 
properties of the RNG introduced in Sec.^ Instead, we seek other kinds of correlations using the random walks test. 
The random walks test proved sensitive and powerful for revealing correlations in RNGs. In particular, correlations in 
the shift register RNG were found ^S^l and explained [3^l40l| using the random walks tests. In addition, the random 
walks test is a useful tool for analysis: if it fails, it gives the opportunity to understand the nature of the correlations 
for a particular RNG [ilf . 

There are several variations of the random walks test in different dimensions ^3 • We consider the one-dimensional 
directed random walk model [s^: a walker starts at some site of an one-dimensional lattice, and at discrete times «, he 
either takes a step in a fixed direction with probability /i or stops with probability 1 — /i. In the latter case, a new walk 
begins. The probability of a walk of length n is P(n) = — /^), and the mean walk length is (n) = 1/(1 — yit). We 
note that the Ising simulations using cluster updates with the Wolff method are closely related to the random walk 
problem Namely, the mean cluster size in the Wolff method equals the mean walk length for /i — tanh( J/fc^T), 
where J is the strength of the spin coupling and T is the temperature. 

Figure 121 shows the correlations in the RNG found by the random walks test. We applied 100 chi-square tests. Each 
test performed n — 10^ random walks with /x = 1/2 for a generator with M = (g 5), w = 32 and s = 1. The result of 
the test with /i = 1/2 is independent of s because only the first bit of the RNG is taken into account. For each test, 
the value 5Pi = (Yj — npi)/ (npi) was calculated for all walk lengths I < 7. Here pi is the theoretical probability of the 
walk length I for uncorrelated random numbers, and Yi is the simulated number of walks with length I. We note that 
correlations can be found only for a large number of random walks (see Table lyi)) . and no correlations are found even 
for n = 6 ■ random walks. 

These correlations can be explained as follows. There are 32 five-bit sequences, and they do not have the same 
frequency of appearing in the RNG output. We consider one of them, for example, 10011. Let X = (0, x (0,1] 

(0) 

and Y = (i, 1] X (0, 1], i.e., X and Y are the left and the right halves of the torus. Let x be the initial point ( "o)) 

of the generator. For the first bits of the first five outputs of the generator to be 10011, it is necessary and sufficient 
to have x S ^looii = Y n R-^{X) n R-^{X) n R-^{Y) n R-'^{Y). Here, R is the action of the cat map. The set 
Zioou consists of filled polygons. Each polygon can be calculated exactly. The area S'(Ziooii) equals the probability 
for the first five outputs of the generator to be 10011. This shows that the nature of the correlations is found in the 
geometric properties of the cat map. 

Figure 13 (the left picture) represents the polygons corresponding to the subsequences of length three for the cat 
map with M = (35). Each set of polygons, e.g., Zqiq = X D i?~^(y) n R~'^{X), represents the region on the torus 
for the first initial point of the RNG and is drawn with its own color. The right picture represents the subsequences 
of length five for the cat map with AI = (J j)- Here, each set of polygons represents the regions on the torus for the 
third point of the generator, e.g., Zqiooi — R^^{X) n R~^(Y) D X O R{X) n R'^{Y), and is drawn with its own color. 
Of course, 5(^01001) = '5'(^oiooi) — P(OIOOI) because the cat maps are area preserving. Therefore, the choice of 
pictures of Zi or pictures of Zt is unimportant if we only want to calculate the areas. Thus, the geometric structures 
in Fig. Oshow the regions of Zqoo, ■ • ■ , ^111 (the left picture) or Zooooo, ■ ■ ■ 1 ^11111 (the right picture) and illustrate 
the geometric approach to calculating the probabilities. 
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FIG. 2: The deviation 5Pi of the probabihty of a walk length / from the value for uncorrelated random numbers versus walk 
length I. The mean and the variance for 5Pi are represented for 100 chi-square random walk simulations. See the text for the 
details. 




FIG. 3: (Color online) Left: the regions on the torus for the first initial point of the RNG described in Sec. ^ with M — (35). 
These regions correspond to the sequences 000, 001, 010, Oil, 100, 101, 110 and 111 of the first bits generated by the RNG. Each 
region is drawn with its own color. Right: the regions on the torus for the third point of the RNG described in Sec. UTI with 
M = {II). These regions correspond to the sequences of length 5 of the first bits generated by the RNG. Each region is drawn 
with its own color. 

The exact areas S'(^ooooo), • ■ • ,5'(^iiiii) can be easily calculated for various toral automorphisms. We prove the 
following geometric propositions: 

1. In any case, every subsequence of length 3, 2 or 1 respectively has the same probability 1/8, 1/4 or 1/2. 

2. If fc = Tr(M) is an odd number, then every subsequence of length 4 has the same probability Pq = 1/16. 

3. If k is even, then the probability of the subsequence 0000 depends only on the trace k of matrix M of the cat 
map. It equals P = Pq ■ k^/{k'^ — 1), where Pq = 1/16. 

The line of reasoning is presented in i45|- Of course, the probability of the subsequence 0000 automatically gives the 
probabilities of all other subsequences of length 4. We note that if k is odd, then ideal (2) is inert (see Appendix I A 3|l . 
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and the inert case is the easiest for exact analysis of the RNG period (see Appendices fXl and IB|) . The probabihties of 
the subsequences of length 5 for maps with odd traces and of the subsequences of length 4 for maps with even traces 
are calculated exactly and shown in Tabled It can be conjectured from Table Ivl that if k is odd, then the probability 
of the subsequence 00000 of length 5 equals Pq ■ {I + l/(3fc2 - 6)), where Pq = 1/32. 

TABLE V: The probabilities of subsequences for different cat maps, characterized by the trace k. 



k 


P(0000)/Po 


k 


P(0000)/Po 


k 


P(00000)/Po 


4 


16/15 


30 


900/899 


3 


22/21 


6 


36/35 


32 


1024/1023 


5 


70/69 


8 


64/63 


34 


1156/1155 


7 


142/141 


10 


100/99 


36 


1296/1295 


9 


238/237 


12 


144/143 


38 


1444/1443 


11 


358/357 


14 


196/195 


40 


1600/1599 


13 


502/501 


16 


256/255 


42 


1764/1763 


15 


670/669 


18 


324/323 


44 


1936/1935 


17 


862/861 


20 


400/399 


46 


2116/2115 


19 


1078/1077 


22 


484/483 


48 


2304/2303 


21 


1318/1317 


24 


576/575 


50 


2500/2499 


23 


1582/1581 


26 


676/675 


52 


2704/2703 


25 


1870/1869 


28 


784/783 


54 


2916/2915 


27 


2182/2181 



The probabilities can thus be approximated as P/Pq = 1 + Bk~^ for large k, where Po — 2~" for subsequences of 
length n = 4,5. Here B = 1 when k is even and n = A; B = 1/3 when k is odd and n = 5. We conclude that the 
deviations found by our implementation of the random walks test will vanish as the trace k increases. 

Table IVTI shows that using rotation in the RNG output (see Sec. Ill D|l results in suppressing correlations found by 
the random walks test. This is not surprising, because even the one-bit random walks test with /i = 1/2 deals with 
the ensemble of cat maps when the rotation is used. 

V. DISCUSSION 

In this paper, we have proposed a scheme for constructing a good RNG. The distinctive features of this approach are 
applying the ensemble of cat maps while taking only a single bit from the point of each cat map and applying methods 
that allow analyzing both the properties of the periodic orbits and the statistical properties of such a generator both 
theoretically and empirically We have seen that the algorithm in Sec.^can generate sequences with very large period 
lengths. Although essential correlations are always present and important statistical deficiencies are found, a good 
algorithm with proper initialization can minimize them. The best generators created by this method have statistical 
properties that are not worse, and speed is slightly slower than that of good modern RNGs. The techniques used 
allow calculating the period lengths and correlation properties for a wide class of sequences based on cat maps. 



TABLE VI: Left: results of the random walks test (i.e. 100 chi-square random walks simulations) for different n and for 
H = 1/2, m = 32 and s = 1. Right: results of the random walks test for different fj, and s and for n = 10^ and m = 32. Here 
s 7^ 1 means using rotation in the RNG output (see Sec. IllDl . Actually, the same one-bit random walks test is used because 
M = 1/2. 



n 


P{K'+) 


PiK'-) 


Result 




s 


P{K'+) 


P{K'~) 


Result 


10* 


0.746106 


0.594428 


PASSED 


1/4 


1 


0.235776 


0.882413 


PASSED 


3- 10* 


0.341899 


0.675728 


PASSED 


3/4 


1 


0.299174 


0.613382 


PASSED 


6- 10* 


0.307433 


0.694282 


PASSED 


1/8 


1 


0.087016 


0.770549 


PASSED 


10^ 


0.018332 


0.966717 


UNCERTAIN 


1/16 


1 


0.67967 


0.920998 


PASSED 


3- lO'^ 


0.0001594 


1 


FAILED 


1/2 


2 


0.605407 


0.344068 


PASSED 


6 ■ 10' 


0.0001378 


1 


FAILED 


1/2 


3 


0.527088 


0.645272 


PASSED 


lO*' 





1 


FAILED 


1/2 


4 


0.558105 


0.360828 


PASSED 
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Future modifications and enhancements are possible, and we currently recommend the generators GM19-SSE, 
GM31-SSE and GRI-SSE for practical use. Program codes for the generators and for the proper initialization can 
be found in and the generator details are discussed in Appendix |E| We would appreciate any comments on user 
experiences. 



VI. ACKNOWLEDGMENTS 



We are grateful to the anonymous referee for the critique and questions that allowed essentially improving the 
content of the paper. This work was supported by the US DOE Office of Science under contract No. W31-109-ENG- 
38 and by the Russian Foundation for Basic Research. 



APPENDIX A: PERIODIC ORBITS OF THE CAT MAPS ON THE 2"x2" LATTICE 



In this section, we review the key arithmetic methods for studying orbit periods that are described in detail in |29l | . 
Some of the results are presented in this appendix in a more general form. The notation is discussed briefly; the 
details and proofs on the formalism of quadratic integers and quadratic ideals can be found in . 



1. The dynamics of the cat map and rings of quadratic integers 

We consider the unit two-dimensional torus (the square (0, 1] x (0, 1] with the opposite sides identified). We take a 
cat map M = " e SL2(L), which acts on a lattice g x g on the torus, where g = 2". The elements of M 

\y 777,21 / 

are integers, detM — 1 and |fc| > 2, where k = Tr(M). 

For any given trace fc > 2, there exists a unique map M G 5*272 (Z) such that the connection between the properties 
of periodic orbits of the automorphism and the arithmetic of quadratic integers is the most natural. Indeed, we 
consider a matrix M such that 

A = mil + T77121, j-^^^l 
Ar = TOi2 + Tm22. 

Here r is the base element of the ring of quadratic integers Rd = {a + br : a,b £ Z} that contains A. This means 
that 377 G Z : fc^ — 4 = n'^D, where _D is a squarefree integer and r = a/D for Z? ^ 1 (mod 4); t — ^(1 + VI)) for 
D = 1 (mod 4). 

It easily follows from ljAl|l that x' + y'r — \{x + yr) is equivalent to (^,) = M(^) for any x,y,x',y'. Indeed, 
X{x + yr) = Xx + (AT)y = (mux + mi2y) + (jn2ix + 77122?/)t = x' + y'r. The action of the map M corresponds 
to multiplication by the quadratic integer A, while the action of M^^ corresponds to multiplication by A^^. Hence, 
we can choose either of the two eigenvalues A = (/c ± Vfc^ — 4)/2, e.g., the largest one, because the exact choice is 
unimportant for studying orbit periods. 

Generally speaking, there are infinitely many maps in 5^2 (^) that have identical eigenvalues, and not all the 
maps are related by a canonical transformation (share the same dynamics). But arguments presented in [29| strongly 
suggest that they still share the same orbit statistics. 



2. Invariant sublattices on the torus and factoring quadratic ideals 

We note that each element of Rjj represents some point of Z^. Let A be a quadratic ideal. We say that ^ = 
T] (mod A) if — rf) € A. We consider the principal quadratic ideal generated by g: (g) = {ag + bgr : a,b G Z}. It 
corresponds to the set of points of a square lattice with the side g. Then the period of an orbit containing the point 
C/P smallest integer T such that z = z (mod (g)). Here x and y are integers, and z = x + yr. 

Each quadratic ideal A is associated with some sublattice of Z^. Because A is a unit, the sublattice is invariant with 
respect to multiplication by A: \A = A. Since we are interested in invariant lattices on the unit two-dimensional torus, 
we consider only those sublattices of Z that are invariant under an arbitrary translation (^^), where a, G Z. These 
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sublattices correspond to quadratic ideals that divide (g). Factoring the ideal (g) thus yields invariant sublattices on 
the torus. 



3. The classification of prime ideals and the orbit periods for the 2"x2" lattice 

We consider 2"x2" lattices on the torus. Because (g) = (2") = (2)", it is sufficient to have the ideal factorization 
of (2). We recall that the ideal (2) is said to be inert if (2) is already a prime ideal; it is said to be split if (2) = P1P2, 
where Pi and P2 are prime ideals; it is said to be ramified if (2) = P^, where Pi is a prime ideal. The ideal (2) is 
inert for _D = 5 (mod 8), split for D = 1 (mod 8) and ramified for D ^ 1 (mod 4). 

It follows that if the trace k is odd, then (2) is inert; if fc = (mod 4), then (2) is ramified. Indeed, for odd k, we 
have fc^ - 4 = 5 (mod 8) ^ D = 5 (mod 8); for fc = (mod 4), we have (fc^ _ 4)/4 = = 3 (mod 4) ^ D = 
3 (mod 4). For k = 2 (mod 4), we obtain (fc^ — 4)/4 = nlD = (mod 4), i.e., all three possibilities (inert, split or 
ramified ideal (2)) can occur. 

Let T„ denote the period of any of the free orbits for g = 2" and denote the period of those ideal orbits for 
g = 2" that do not belong to the sublattice f x|. We recall that an orbit belonging to a given lattice Z /gZ is 
called an ideal orbit if it belongs to some ideal A such that A\{g) and Aj^ (1). Otherwise, it is called a free orbit. 

The behavior of periodic orbits on the 2x2 lattice follows from Propositions B1-B3 in |29|. Namely, we have the 
following: 

• If (2) is inert, then either Ti = 3 or Ti = 1; all orbits are free. 

• If (2) is split, then Ti = T{ = 1; there are two ideal orbits and one free orbit. 

• If (2) is ramified, then Ti = 2 and = 1; there is an ideal orbit and a free orbit (it is also possible that Ti = 1 
and T[ — 1] there are two free orbits and an ideal orbit). 

To determine the structure of periodic orbits on the 2"x2" lattice, we prove the following theorem. 
Theorem. 

1. For all ri, either r„+i = 2r„ or T„_|_i = T„. 

2. For aU n, either T'^ = T„ or T.;^ = T^-i. 

3. For aU n > 3, T„ ^ T„_i =^ T„+i ^ T„. 

4. If n > 4, Tn Tn-i, and = r„/a, where a G {1, 2}, then T^_|_i — Tn+i/a. 

This theorem generalizes Propositions CI and C2 in (291]. The line of reasoning is presented in Appendix IdI 

Therefore, knowing T„ and for small n suffices for determining the orbit statistics for all n. There always exist 
ni, n2 and ng such that r„ = ri2"-"i and T/, = ri2"-"2 for all n > rig. 

If (2) is inert, then every ideal that divides (g) has the form (2)''. Therefore, each ideal orbit belongs to the 
2"^^x2"^^ sublattice and coincides with a free orbit for some sublattice 2''x2'', where r < n. We now find the 
number of free orbits in the inert case. There are 2^"— 1 points on a lattice. The ideal orbits contain 2^"^^ — 1 points. 
Consequently, there are (2^" - 22"-2)/r„ = 3 • 22"-^7T„ free orbits. 

We suppose that the typical inert case occurs, i.e., T„ = 3 • 2"^^. Then the phase space is divided into the following 
regions: 

• 3/4 of the phase space is swept by 2" trajectories of period T„, 

• 3/16 of the phase space is swept by 2^~-^ trajectories of period T„_i — T„/2, 

• 3/64 of the phase space is swept by 2"^^ trajectories of period T„_2 = Tri/4, 

• and so on. 

All such statements hold as long as the trajectory length exceeds just a few points. Therefore, on one hand, cat map 
orbits have huge periods; on the other hand, the number of orbits is sufficiently large (see Theorem 2 in Appendix IB)| . 
Both these properties are important for our construction of the RNG. 
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APPENDIX B: THE RNG PERIOD 

In this section, we find the periods of the generators in Sec. |n] and Sec. IIIDI As a result of Appendix ^ and 
Appendix^ the RNG period can be obtained for arbitrary parameters of the map and lattice. 

Theorem 1. If g = 2™, then the period T of the sequence {a^"-'} in Sec. UTI equals the period T„j of free orbits of 
the cat map for the overwhelming majority of RNG initial conditions. 

Proof. 

1. At least one of the initial points ( (o)) belongs to a free orbit. Indeed, the probability of this in the inert case 
equals (1 - 4""*). 

2. Therefore, T is not less than r,„. Indeed, T is not less than the period of the sequence of first bits of xf^ , xf'^ , . . . 
for each i. But the period of the sequence of first bits of points of the cat map orbit is equal to the orbit period 
for the vast majority of orbits. The probability of the opposite is tiny provided that the orbit is not too short. 

3. Finally, T is not larger than Tm- Indeed, the period of each cat map orbit divides T^. 

Example. In the typical example of the inert case, where M = (37), we obtain = 3 • 2™~^. This fact was 
also tested numerically as follows. First, the initial conditions were set randomly. Second, the period of {a^")} was 
accurately found numerically. This operation was repeated 1000 times for m = s = 14. Each time the period of {a'^*-'} 
turned out to be 6144 = 3 • 2^^ . To check the period numerically, we first check whether the whole state of the RNG 
(not only the output) coincides at the moments and T and then verify that a smaller period (which could possibly 
divide T) does not exist. 

Theorem 2. The probability that two arbitrary points of the 2™x2'" lattice on the torus belong to the same orbit 
of the cat map equals 9/(7 • 2™+^). The probability that s arbitrary points of the lattice do not belong to s different 
orbits of the cat map (i.e., two of the points belong to the same orbit) is 9s(s— 1)/(7 • 2™+^). 

The proof of Theorem 2 is straightforward. Of course, both these probabilities are tiny if m is sufficiently large. 

Theorem 3. If g = 2™ and s|T„j, then the period T of the sequence {fe^^"^} in Sec. Ill Dl equals r,„ for the 
overwhelming majority of RNG initial conditions. 

Proof. 

1. Because s\Tm, we have &i+T,n = &i for all i. Therefore, T\Tm- 

(0) 2.(0) 

2. If s does not divide T, then Vi G {0, 1, . . . , s — l}3j G {0, 1, . . . , s — 1}, j 7^ i, such that ( joi) and ( L-.) belong 

Vi y) 

to the same orbit of the cat map. It follows from Theorem 2 that this event is highly improbable. Therefore, 
s\T. 

3. Because T is a period of {fc*^"^} and s\T, we have &i+T — &i ^ o-i+T ~ ai for all i. Therefore, Tm\T. 

The above theorems show that the period calculations for the sequences {a'^"^} and {&^")} are reliable in the general 
case, because the chance of the period dependence on the initial state is exponentially small. But it is a desirable 
property that the period does not depend on any conditions at all. The proper initialization (see Appendix IE|) 
guarantees that (i) at least one of the initial points belongs to a free orbit of the cat map; and (ii) no pair of initial 
points belongs to the same orbit of the cat map. Therefore, both the periods of {a*^"-'} and of {6'^"^} are guaranteed 
to equal Tm provided the initialization in Appendix ^ is applied. 

The above theorems and considerations hold for g = 2™. In the other case, when g = p is a, prime, it follows from 
finite field theory that the period of any orbit of the matrix transformation is equal to — 1 provided the polynomial 
fix) = x^ — kx-\-q is primitive modulo p. The methods for good parameter and initialization choice for such generators 
are also presented in Appendix IeI for the generators GM19 and GM31. A similar argument as in Theorems 1 and 3 
shows that in this case (i) the period of the sequence {a^"^} equals — and (ii) the period of the sequence {6^"^} 
is divisible by p^ — 1, i.e., rotation cannot decrease the period of such a generator. 

APPENDIX C: ORBITS, NORM AND CORRELATIONS BETWEEN ORBITS 

In this section, (i) we show that the norm modulo g is the characteristic of the whole orbit; (ii) we find the number 
of orbits of each norm modulo g and discuss how symmetries affect the norm; and (iii) we find the linear congruential 
dependences between orbits. The consideration holds for maps with q = 1 on a 2™x2'" lattice. 
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1. Orbits and norm 

We recall that the norm of a quadratic integer a = a + b\fD is simply an integer N{a) — aa* = — b^D. If (2) 
is inert, then the quadratic integer x + yr, where r = ^-^j^, represents the point (p, i.e., iV(y) = + a;?; — ^^j-^y^- 
A cat map preserves the value of N{^) (mod 2™), because the action of a cat map can be described as x' + y'r — 
X{x + yT) ( mod (2™)), where A is a matrix eigenvalue and N{X) = 1. Therefore, the norm modulo g is a characteristic 
of the whole orbit. We note that for a point on a free orbit, either a; or y is odd, consequently the norm is also an 
odd number. 

We prove that if the period of free orbits is T = 3 • 2™^^, then for each N — 1,3,..., 2™— 1, there are exactly two 
orbits that have the norm N (these two orbits are symmetrical, i.e. the second one contains the points (2mZy")j where 
(^") are the points of the first one). Indeed, there are exactly 2™ free orbits that occupy T • 2™ — 2^™ — 22™-2 points. 
On the other hand, there is a method for obtaining two symmetrical orbits having any odd norm. We note that other 
possible symmetries (e.g., symmetries considered in preserve the norm modulo 256. Moreover, in most cases, 
they preserve the norm modulo 2™^^ or modulo 2™^^. 

2. Correlations between orbits 

We consider a pair of free orbits with the norms A^i and A^2- The set A — {1, 3, . . . , 2™— 1} is a group under 
multiplication (it is called the modulo multiplication group); hence, there exists t G A such that A^i = tN2 (mod 2™). 
It is known that for the equation = t (mod 2™) to have a solution k G A, it is necessary and sufficient to have 
t = 1 (mod 8). Thus, if A^i = N2 (mod 8), there exists k such that iVi = k^N2 (mod 2™). If (^'^) are the points of 

the orbit of norm then (^^") (mod (2™)) are the points of the orbit of norm iVi = k^N2 (mod 2™), in the same 
order. But there may be large shift between the values of different orbits. 

Thus, the case iVi = N2 (mod 8) is dangerous, because there may be correlations between orbits. The points of 
the first orbit are connected to the points of the second orbit with a linear congruential dependence. The parameter 
k may be interpreted as a random odd number. 

APPENDIX D: PROOF OF THE THEOREM IN APPENDIX FOl 

Proposition 1. T„ is the least integer such that A^" = 1 (mod (2")). In particular, any free orbit has the same 
period. 

Proof. We suppose that the period of an orbit containing the point z is T. Then X^z = z (mod (2")) 
z{X^-\)g (2"). If the orbit is free, then z ^ P for any ideal P such that P| (2" >, P 7^ (1). Therefore, (A^-l)e (2"). 
Proposition 2. Tn\Tnj^\. 

Proof. Indeed, A^"+i = 1 (mod (2"+^)) A^"+i ee 1 (mod (2")) => r„+i = mT„, where m e N. 
Proposition 3. For all n, either T„+i — 2Tn or T„+i = r„. 

Proof. Because (A"^" - 1) S (2"), we have (A^" + 1) = (A"^" - 1) + 2 e (2). Consequently, A^^" - 1 = 
(A^- - 1)(A^" + 1) e (2"+i), i.e., either P„+i = 2r„ or T„+i = T„. 
Proposition 4. If n > 3 and T„ ^ T„_i, then T^+i ^ r„. 

\ ~ ]\ ,'' =^ A"^"-! = 1 + z- 2"-i, where z ^ (2). 

X "-^ ^ 1 (mod (2")) 

Squaring the last equation, we obtain A^^"-i = 1 + z • 2" + • 22"-^ ee 1 + z • 2" (mod (2"+^)) for n > 3. Hence, 
A2T„-i ^ 1 (j^od (2"+i)) ^ r„+i ^ r„. 

Proposition 5. If (2) is split, then for all n, T!^ — T„. In particular, T!^ is the same for all ideal orbits that do not 
belong to the sublattice 2"~^ x 2"~^, no matter what the ideal is. 

Proof. Because (2) is split, we have (2) = PiP2- Let T and S be the smallest integers such that A"'" = 1 ( mod P") 
and A"^ = 1 (mod PJ*). We prove that T = S*. First, we note that Pi and P2 are conjugate ideals, i.e.. Pi = Pj*. 
We assume T = S + R and P > 0. Taking the conjugate of the congruence A'^ = 1 (mod P^^), we obtain A*"^ = 
1 (mod Pf), where A* = A^\ Therefore, X^ X*^ X^ = = 1 (mod Pf) ^ A^ = 1 (mod Pf), i.e., there exists an 
integer I > such that R = IT. Because T = S + lT,we have I = 0^T = S. 

Let z belong to an ideal orbit of length P,' and z G P2 , z ^ P2^^, where k e {1, 2, . . . , n}. Then T,' and P„ are 
the smallest integers such that A"^" = 1 (mod P^P^~'') and A"^" = 1 (mod Pi"P2"). Therefore, r;;|T„. On the other 
hand, A^" = 1 (mod Pf) ^ A^- = 1 (mod P2") ^ A^- = 1 (mod PfP^f ), i.e., r„|r;;. Therefore, = r„. 
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Proposition 6. If (2) is ramified, then for all n, either = Tn or = T„_i. 

Proof. We have (2) = . We consider an orbit belonging to P. We now show that the orbit period is either T„ 

or r„_i. 

' A"^"-! = 1 (mod (2"-i)), 
< A^" = 1 (mod (2"-i)P), ^ r„_i|T^|T„. 
^ A^" = 1 (mod (2")). 

Using Proposition 3, we complete the proof. 

Proposition 7. If (2) is ramified, n > 3 and r„ = 2r„_i, then r,;+i = 21;',. 
Proof. Let T ~ T„_i. Then we have 

f A^ = 1 (mod A), 

\ A"^ ^ 1 (mod AP), ^ ' 

where A ^ (2"-i> for T!^ ^ r„ and A = (2"-i)P for T/, ^ T^-i. In any case, (2"-i)|A, AP|(2">. It follows from jDll) 
that A^ 1 + z, where z e A, z ^ AP. Hence, A^^ = 1 + 2z + z^. We note that 2z G ((2)A), 2z ^ ((2)AP), and 
^2 g ^2"+i) ^ G ((2)AP) for n > 3. Therefore, 

A2^ ^ 1 (mod (2)A), 

A2^ ^ 1 (mod {2)AP). ^ ' 

If = T„, this means that T/j_|.i ^ T,i 7,'+! = Pn+i- In the case where T/^ = T„_i, we have T/j+i = Tn- In any 
case, T/j_|_]^ = ^T^^. 



APPENDIX E: REALIZATIONS AND ALGORITHMS 



1. RNG realizations in C language and in inline assembler, speed of realizations 

In this section, we present efficient algorithms for several versions of the RNG introduced in Sec^l I^i particular, 
GS (cat map Generator, Simple version), GR (cat map Generator, with Rotation), GRI (cat map Generator, with 
Rotation, with Increased trace), GM (cat map Generator, Modified version). The parameters and characteristics for 
these generators can be found in Tabic IVTll and the results of stringent statistical tests in Sec. IIIII For comparison, 
both in Table Ivnl and in Sec. lIIIBl w e also test the standard UNIX generators randO, rand48() and random () and 
the modern generators MT19937 [li|. MRG32k3a 15] and LFSR113 [13 (see Sec. UlTBl for details on them). 

Most of our generators are speeded up using Eq. (O instead of Eq. Also, the Streaming SIMD Extensions 2 
(SSE2) technology, introduced in Intel Pentium 4 processors allows using 128-bit XMM-registers to accelerate 
computations. A similar technique was previously used for other generators The SSE2 algorithms for our 

generator are able to increase performance up to 23 times as compared with usual algorithms (see Table IVIl|l . 

The algorithms for GS, GRI and GM19 are shown in Table IVlTll Table Hxl illustrates the key ideas for speeding up 
cat-map algorithms using SSE2. We use the GCC inline assembler syntax for the SSE2 algorithms. The action of the 
fast SSE2 algorithms shown in the left column are equivalent to the action of the slow algorithms shown in the right 
column. 

The complete realizations for all RNGs can be found in GM31-SSE is the only algorithm here that exploits 
64-bit SSE-arithmetic for calculating Eq. jS)). We must also note that the algorithms that exploit the SSE2 command 
set work properly for Pentium processors starting from Pentium IV. Therefore, some of our codes are not immediately 
portable. Even the AMD's implementation of SSE2 is based on a slightly different command set. 



2. Initialization of generators 



The proper initialization is very important for a good generator. 

For the generators GS, GS-SSE, GR-SSE, GRI, GSI-SSE and GRTSSE we use the following initialization method: 

^(0) 

• Norms of all points should be different modulo 256. In particular, this guarantees that the initial points ( Jo)), 

i = 0, 1, . . . , (s — 1) belong to different orbits of the cat map, and that none of the symmetries may convert one 
orbit to another (see Appendix O. 
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TABLE VII: Characteristics and parameters for several versions of the RNG based on the ensemble of cat maps (see Sec. ^ 
and for other generators (last six entries). Here "CPU-time" means the CPU time (in seconds) needed to generate 10* uniform 
random numbers on a 3.0 GHz Pentium 4 PC running Linux. This parameter characterizes the speed of the generator. 
Generators may be used only when the application needs not more that T random numbers, where T is the RNG period. 



Generator 


g 


s 


k 


q 


Rotation 


SSE2 


Period 


CPU-time 


GS 


232 




32 


3 


1 






3.2 • lO'' 


55.4 


GS-SSE 


232 




32 


3 


1 




+ 


3.2 • lO'' 


2.49 


GR-SSE 


232 




32 


3 


1 


+ 


+ 


3.2 • lO'* 


2.79 


GSLSSE 


232 




32 


11 


1 




+ 


3.2 • lO'* 


3.66 


GRI 


232 




32 


11 


1 


+ 


— 


3.2 • lO'* 


78.2 


GRLSSE 


232 




32 


11 


1 


+ 


+ 


3.2 • lO'* 


4.03 


GM19 


219 _ 


1 


32 


6 


3 


+ 


— 


2.7 ■ 10" 


120.5 


GM19-SSE 


2i9 _ 


1 


32 


6 


3 


+ 


+ 


2.7 ■ 10" 


6.11 


GM31-SSE 


231 _ 


1 


32 


7 


11 


+ 


+ 


4.6 • 10^* 


8.86 


RAND 
















2.1 • 10" 


2.48 


RAND48 
















2.8 • 10" 


4.64 


RANDOM 
















3.4 ■ 10^° 


1.88 


MT19937 
















4.3 ■ 10*°^ 


2.45 


MRG32k3a 
















3.1 • 10" 


11.14 


LFSR113 
















1.0 ■ lO^"* 


2.98 



• At least one point should belong to a free orbit, i.e., at least one of the coordinates x or y should be an odd 
number. This guarantees that the period length is not smaller than (see Appendix IB|| . 

We choose the parameters k and q for the generators GM19 and GM31 such that the polynomial f{x) = x'^—kx+q 
is primitive modulo p, where p — 2^^— 1 for GM19 and p = 2^^ — l for GM31. Therefore, the actual period of the 
generator is p^—1. 

To construct the initialization method for GM19 and GM31, we use the "jumping ahead" property, the possibility 
to skip over terms of the generator. In other words, we utilize an easy algorithm to calculate Xn quickly from Xq and 
Xi, for any large n. We choose the following initial conditions: xf''^ = x^j^"^ — Xi^A+i^i = 0, 1, . . . ,31. Here we 

follow the notation in Sec. ITlland A is a value of the order of (p^ — 1)/32. We recommend to choose A randomly; at 
least A should not be chosen very close to the divisor of — 1 or to a large power of 2. We recommend using less than 
A random numbers in applications that use GM19 and GM31. The values of A are approximately 32 times smaller 
than the periods in Table IVIll 

The initialization routines for all generators can also be found in l43| . 



[1] K.S.D. Beach, P.A. Lee, P. Monthoux, Phys. Rev. Lett., 92 (2004) 026401. 

[2] D.P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University Press, 
Cambridge, 2000). 

[3] S.C. Pieper and R.B. Wiring, Ann. Rev. Nucl. Part. Sci., 51 (2001) 53. 
[4] A. Liichow, Ann. Rev. Phys. Chem., 51 (2000) 501. 
[5] A.R. Bizzarri, J. Phys.: Cond. Mat., 16 (2004) R83. 

[6] D.E. Knuth, The art of Computer Programming, Vol. 2, (Addison- Wesley, Cambridge, 1981). 
[7] P. L'Ecuyer, Ann. of Oper. Res., 53 (1994), 77. 

[8] R.P. Brent and P. Zimmermann, in Lecture Notes in Computer Science, Comp. Sc. and its Appl. - ICCSA 2003 (Springer- 
Verlag, Berlin, 2003), 1. 

[9] R.R. Coveyou and R.D. MacPherson, J. ACM 14 (1967) 100; G. Marsaglia, Proc. Nat. Acad. Sci. USA 61 (1968) 25. 
[10] S.W. Golomb, Shift Register Sequences, (Holden-Day, San Francisco, 1967). 
[11] A.M. Ferrenberg, D.P. Landau, Y.J. Wong, Phys. Rev. Lett., 69 (1992) 3382. 
[12] P. Grassberger, Phys. Lett., 181 (1993) 43. 
[13] F. Schmid, N.B. Wilding, Int.J.Mod.Phys., C 6 (1995) 781. 

[14] M. Matsumoto and T. Nishimura, ACM Trans, on Mod. and Comp. Sim., 8 (1998) 3. 

[15] P. L'Ecuyer, Oper. Res., 47 (1999) 159. 

[16] P. L'Ecuyer, Math, of Comp., 65 (1996) 203. 



17 



TABLE VIII: Codes in ANSI C language for the generators GS, GRI and GM19. 



const unsigned long halfg=2147483648; 
unsigned long x[32],y[32]; char rotate; 



//- 



Generator GS 



unsigned long GS(){ 

unsigned long i , output=0 ,bit=l ; 
for(i=0;i<32;i++){ 
x[i]=x[i]+y[i] ; 
y[i]=x[i]+y[i] ; 

} 

for(i=0;i<32;i++){ 

output+=((x[i]<halfg)?0:bit) ; bit*=2;} 
return output ; 

} 



II- 



Generator GRI 



unsigned long GRI()-[ 

unsigned long i , oldx, oldy , output=0 ,bit=l ; 
oldx=x[31]; oldy=y[31]; 
for(i=31;i>0;i— ){ 

X [i] =4*x [i-1] +9*y [i-1] ; 

y [i] =3*x [i-1] +7*y [i-1] ; 

}; 

X [0] =4*oldx+9*oldy ; y [0] =3*oldx+7*oldy ; 
for(i=0;i<32;i++){ 

output+=((x[i]<halfg)?0:bit) ; bit*=2;} 
rotate++; return output; 



II- 



Generator GM19 



const unsigned long k=14; 

const unsigned long q=15; 

const unsigned long g=524287; 

const unsigned long qg=7864305; 

const unsigned long half g=262143; 

unsigned long x [2] [32] ; char new, rotate; 

unsigned long GM19(){ 

unsigned long i , output=0 ,bit=l ; 
char old=l-new; 
for(i=0;i<32;i++) 

x[old] [i] = (qg+k*x[new] [i]-q*x[old] [i])"/.g; 
for(i=0;i<32;i++){ 

output+=((state->x[old] [(256+i-rotate)y.32] <half g) ?0 :bit) ; 

bit*=2; 

> 

new=old; rotate++; return output; 
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TABLE IX: Equivalent realizations for several algorithms with inline assembler code for Pentium 4 processor (left column) and 
ANSI C language (right column). First row presents the main part of the GRI algorithm. Second row presents the packing 16 
high bits of 16 integers into one integer. These or similar equivalences are used in constructing the SSE2 algorithms for any of 
the discussed RNGs 0. 



unsigned long x[4],y[4]; 
[ ] 

asmC'movaps (7,0) //.y.xmmOXn" \ 
"movaps C/.l) ,'/.y,xmml\n" \ 
"paddd •/.•/.xmml,y;/.xminO\n" \ 
"paddd y.y.xmml,y.y.xmmO\n" \ 
"movaps y,y,xmm0,y,y,xmm2\n" \ 
"pslld $2,y.y.xmin0\n" \ 
"paddd y.y.xmml,y.y.xminO\n" \ 
"movaps y.y.xmmO, (y.O)\n" \ 
"psubd y,y,xmm2,y,y,xmm0\n" \ 
"movaps y.y.xmmO , (y.l) \n" \ 
"": :"r"(x),"r"(y)); 


unsigned long i ,newx [4] ,x [4] ,y [4] ; 
[ ] 

for(i=0;i<4;i++){ 

newx [i] =4*x [i] +9*y [i] ; 
y [i]=3*x[i]+7*y [i] ; 
x [i] =newx [i] ; 

} 


unsigned long x [16] , output ; 
[ ] 

asm("movaps (y.l) ,y,y,xmmO\n" \ 
"movaps 16(y,l) ,y,y.xmml\n" \ 
"movaps 32(y.l) ,y.y.xmm2\n" \ 
"movaps 48(y.l) ,y.y.xmm3\n" \ 
"psrld $31,y.y.xmm0\n" \ 
"psrld $31,y.y.xmml\n" \ 
"psrld $31,y.y.xmm2\n" \ 
"psrld $31,y.y.xmm3\n" \ 
"packssdw y.y.xmml //.y.xmmOXn" \ 
"packssdw y.y,xmm3,y,y.xmm2\n" \ 
"packsswb y,y,xmm2,y,y,xmm0\n" \ 
"psllw $7, y.y.xmmO \n" \ 
"pmovmskb y.y.xmmO, y,0\n" \ 
"" : "=r" (output) : "r"(x)) ; 


const unsigned long half g=2147483648 ; 
unsigned long x[16] ,i,output=0,bit=l; 

[ ] 

for(i=0;i<16;i++){ 

output+= ( (x [i] <half g) ?0 : bit ; 
bit*=2; 

> 
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